#include "SparseBlas.h"

using namespace SG::Algebra;

//=============================================================================
// SparseBlas

SparseBlas::SparseBlas ()
{

}

SparseBlas::~SparseBlas ()
{

}

#if defined(SG_USE_MPI)
Info_t SparseBlas::setValues (ComplexSpMat_t& mat, const Int_t& i_start, const Int_t& j_start, const SpMat_t& sm)
{
    typedef Triplet<Complex_t> Triplet;

    int exit_status = 0;

    // fetch non-zero elements from the sub-matrix
    std::vector<Triplet> triplets;

    auto M = mat.rows ();
    auto N = mat.cols ();

    auto row_start = mat.row_start ();
    auto row_end   = mat.row_end ();

    auto nnz = mat.nonZeros ();
    // Requests that the vector capacity be at least enough to contain nnz elements.
    triplets.reserve (nnz);

    std::vector<Index_t> outerIndex;
    std::vector<Index_t> innerIndex;
    std::vector<Real_t>  values;

    sm.get_csr (outerIndex, innerIndex, values);

    for (auto i = row_start; i < row_end; ++i)
    {
        auto index_start = outerIndex[i];
        auto index_end   = outerIndex[i + 1];

        for (auto index = index_start; index < index_end; ++index)
        {
            auto j = innerIndex[index];
            auto v = values[index];

            // skip out the out-bounded element
            if (j > N - 1)
            {
                continue;
            }

            triplets.push_back (Triplet (i, j, v));
        }
    }

    // insert non-zero elements into the blocked matrix
    // mat.setFromTriplets (triplets.begin (), triplets.end ());

    ComplexSpMat_t mat2 (M, N);
    mat2.setFromTriplets (triplets.begin (), triplets.end ());

    mat = mat + mat2;

    return ( Info_t ) exit_status;
}

Info_t SparseBlas::dot (const Vec_t& X, const Vec_t& Y, Real_t& R)
{
    int exit_status = 0;

    R = X * Y;

    return ( Info_t ) exit_status;
}

Info_t dot (const ComplexVec_t& X, const ComplexVec_t& Y, Complex_t& R)
{
    int exit_status = 0;

    R = X * Y;

    return ( Info_t ) exit_status;
}
#else
Info_t SparseBlas::setValues (ComplexSpMat_t& mat, const Int_t& i_start, const Int_t& j_start, const SpMat_t& sm)
{
    typedef Eigen::Triplet<Complex_t> Triplet;

    int exit_status = 0;

    // fetch non-zero elements from the sub-matrix
    std::vector<Triplet> triplets;

    auto M = mat.rows ();
    auto N = mat.cols ();

    auto m   = sm.rows ();
    auto n   = sm.cols ();
    auto nnz = sm.nonZeros ();
    // Requests that the vector capacity be at least enough to contain nnz elements.
    triplets.reserve (nnz);

    for (auto j = 0; j < n; ++j)
    {
        // skip out the out-bounded element
        if (j > N - 1)
        {
            continue;
        }

        auto index_start = sm.outerIndexPtr ()[j];
        auto index_end   = sm.outerIndexPtr ()[j + 1];

        for (auto index = index_start; index < index_end; ++index)
        {
            auto i = sm.innerIndexPtr ()[index];
            auto v = *(sm.valuePtr () + index);

            // skip out the out-bounded element
            if (i > N - 1)
            {
                continue;
            }

            triplets.push_back (Triplet (i + i_start, j + j_start, v));
        }
    }

    // insert non-zero elements into the blocked matrix
    // mat.setFromTriplets (triplets.begin (), triplets.end ());

    ComplexSpMat_t mat2 (M, N);
    mat2.setFromTriplets (triplets.begin (), triplets.end ());

    mat = mat + mat2;

    return ( Info_t ) exit_status;
}

Info_t SparseBlas::dot (const Vec_t& X, const Vec_t& Y, Real_t& R)
{
    int exit_status = 0;

    R = X.dot(Y);

    return ( Info_t ) exit_status;
}

Info_t dot (const ComplexVec_t& X, const ComplexVec_t& Y, Complex_t& R)
{
    int exit_status = 0;

    R = X.dot (Y);

    return ( Info_t ) exit_status;
}
#endif 

Info_t SparseBlas::multiply (const SpMat_t& A, Vec_t& X, Vec_t& Y)
{
    Y = A * X;

    return Info_t (0);
}

Info_t SparseBlas::multiply (const ComplexSpMat_t& A, ComplexVec_t& X, ComplexVec_t& Y)
{
    Y = A * X;

    return Info_t (0);
}
